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1.0  Introduction 


In  many  applications,  such  as  the  phased  array  radar  problem,  we  have  an 
observed  J  channel  data  vector  z(n)  at  time  n  which  may  contain  a  desired  signal 
s(n)  and  an  unwanted  additive  disturbance  y.(n)  which  may  contain  additive 
temporally  correlated  noise  (clutter)  c(n)  and  uncorrelated  thermal  white  noise 
w(n).  For  the  detection  problem,  a  decision  must  be  made  between  the  two 
hypotheses 

Ho:z(n)  =  y(n)  n=l,2,...,N 

z(n)  =  s(n)  +  y(n)  n=l,2,...,N.  (1) 


For  high  resolution  radar,  the  spatial  clutter  statistics  deviate  from  the 
Gaussian  model  [1,2].  The  non-Gaussian  nature  of  the  clutter  can  cause  severe 
performance  reduction  in  applications  such  as  space-time  processing  and  bearing 
estimation  using  phased  array  systems.  In  these  applications,  a  covariance  matrix 
is  used  to  estimate  the  unwanted  interference  noise  in  the  data  cell  under  test. 
Ensemble  averaging  is  replaced  by  averaging  over  available  adjacent  secondary 
data  cells  (i.e.,  range  cells  assumed  to  be  free  of  the  desired  signal)  assuming 
independence  between  range  cells.  This  estimator  is  the  sample  matrix  estimator 
and  we  refer  to  this  case  as  that  of  a  spatially-averaged  estimate.  The  presence  of 
the  spatially  non-Gaussian  clutter  causes  a  considerable  increase  in  the  variance  of 
the  sample  covariance  matrix  estimator  as  compared  to  that  for  the  Gaussian  case. 

In  other  applications  such  as  the  airborne  radiometer  [3],  a  time-averaged 
estimate  is  used  by  averaging  over  a  time  sequence  of  data.  In  this  paper,  the 
variance  of  the  spatially-averaged  and  time-averaged  covariance  estimators  are 
described  analytically  in  terms  of  parameters  of  non-Gaussian  models  pertinent  to 
these  applications. 


Two  families  of  distributions  for  SIRP's  have  received  considerable 
attention;  namely,  the  K-distribution  [1] 


fta+lrOt 


r>0 


P,oo0 
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and  the  WeibuU  [2], 

fj^(r)  =  pai"'^exp(-pi^)  r  >  0  P,cx>0  (2) 

where  r(*)  is  the  Eulerian  function  and  K(*)  is  the  modified  Bessel  function  of  the 
second  kind  of  order  a.  The  parameters  P  and  a  are  the  scale  and  shape 
parameters,  respectively.  The  former  is  related  to  the  variance  of  the  quadrature 
components  of  the  clutter,  while  the  latter  controls  the  tails  of  the  distribution. 

The  complex  envelope  y(n,k)  at  time  n  and  range  ceU  k  of  the  compound- 
Gaussian  clutter  model  [2,4]  has  been  considered;  i.e., 

y(n,k)  =  v(n,k)x(n,k)  (3) 

where  the  complex  Gaussian  process  x(n,k)  (the  speckle  component)  modulates  a 
nonnegative  process  v(n,k)  containing  discretes  (texture)  which  is  statistically 
independent  of  x(n4c).  When  v(n,k)  has  a  long  coherence  time,  it  can  be  considered 
to  be  a  random  variable  with  respect  to  k,  but  constant  over  time  n  so  that  eq(3) 
reduces  to  the  case  of  a  spherically  invariant  random  process  (SIRP)  [5,6] 

y(n,k)  =  v(k)x(n,k).  (4) 

Eq(4)  is  known  as  the  representation  theorem  [5]  and  the  probability  density 
function  (pdf)  ^(v)  is  called  the  characteristic  pdf.  It  controls  the  pdf  of  y(n,k). 

The  WeibuU  and  K-distributed  processes  are  both  special  cases  of  the  SIRP's  [6]. 
We  note  that  conditioned  on  v(k),  the  time  sequence  y(n)  is  Gaussian.  This  implies 
that  a  sequence  of  data  from  a  specific  range  cell  is  locally  Gaussian,  but  non- 
Gaussian  when  observed  across  cells.  Furthermore,  y(n,k)  is  non-ergodic;  i.e., 
averaging  over  range  ceUs  where  v(k)  is  a  random  variable  from  ceU-to-cell  does 
not  provide  the  same  statistic  when  averaging  over  time  with  v(k)  constant. 

In  this  paper,  the  variance  of  both  the  spatially-averaged  and  time-averaged 
covariance  estimates  is  derived.  The  special  cases  of  WeibuU  and  K-distributed 
processes  are  noted  specificaUy.  Another  important  contribution  of  this  paper  is 
the  use  of  the  representation  theorem  of  eq(4)  in  the  derivation.  This  theorem 
enables  a  significant  simplification  of  the  procedure. 
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2.0  Derivation 


2.1  Spatially  Averaged  Estimator 

We  first  consider  the  spatially-averaged  covariance  estimator;  i.e.  averaging 
performed  over  range  cells  as  described  above.  Consider  the  J  channel  vector 
y.(n,k)  with  JxJ  covariance  matrix  R=E[y.(n,k)y^^(n,k)]  where  H  denotes  the 
Hermitian  transpose.  Suppressing  the  time  dependence,  let  x(k)  be  a  J  channel 
vector  of  observation  processes  from  the  k***  sample  (or  secondary  range  cell). 
Consider  the  JxJ  covariance  matrix  estimate  defined  as 

^  X(k)y^(k).  (5) 

k=l 

Assuming  statistical  independence  over  k,  ^  is  the  sample  covariance  matrix, 
and  provides  a  maximum  likelihood  estimate  for  Gaussian  processes  [7].  The  ij^ 
element  ^ y  of  ^  is  defined  as 

^j  =  KSyi(k)yJ(k).  (6) 

k=l 

The  estimator  is  unbiased  since  E[f  jj]  =  r^.  Thus  the  variance  and  error  variance 
are  identical.  The  variance  Vj.^  of  f  jj  is  defined  as, 

V,3  =  E{[fij-E(fij)][i^y-E(^ij)r} 

=  -  ECi^ipEC^J)  (7) 

where  the  subscript  s  on  V  denotes  the  spatially  -averaged  case.  Using  (6)  in  (7) 
and  interchanging  the  summation  and  expectation  operations, 

Vrs  =  E  E{yi(m)yJ(m)y*(p)yj(p)} 

^  m=lp=l 

1  K  ^  ^ 

-  E  E  E[yi(m)yj  (m)]E[yi  (p)yj(p)]  (8) 

^  m=lp=l 
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where  m  and  p  denote  the  m*  and  range  cells,  respectively.  Using  the 
representation  theorem  from  eq(4)  and  noting  the  statistical  independence  of  v  and 
Xj,  the  first  expectation  in  eq(8)  becomes 

E  { yi(m)y* (m)y* (p)yj(p) }  =  E[VmVp]  E  { Xi(m)Xj  (m)xi  (p)Xj(p) } .  (9) 

For  zero  mean,  jointly  Gaussian  scalar  random  variables  Xj(*),  Xj(*),  eq(9) 
reduces  to 


E  {  yi(m)y* (m)y* (p)yj(p) }  =  E[VmVp]  (  E[Xi(m)Xj  (m)]E[xi  (p)Xj(p)] 

+  E[Xi(m)x*(p)]E[x*(m)Xj(p)] } 

+  E[xj(m)Xj(p)]E[Xj  (m)Xi  (p)].  (10) 

Furthermore,  when  these  processes  are  complex,  circularly  symmetric,  the 
expectations  in  the  last  term  of  eq(lO)  are  zero,  so  that 

E  {  yi(m)y* (m)y* (p)yj(p) }  =  E[VmVp]  {  E[Xi(m)Xj  (m)]E[Xi  (p)Xj(p)] 

+  E[Xj(m)x*(p)]E[Xj  (m)Xj(p)] } .  (11) 

Assuming  spatial  stationarity  over  range  cells,  eq(ll)  becomes 

E  {  yi(m)yj (m)yi  (p)yj(p) }  =  E[VmVp]  { lr”(0)|2  +  ry(m-p)[rj(m-p)]  *  }  (12) 

where  the  superscript  x  denotes  that  the  correlation  functions  describe  the  Gaussian 
processes  defined  in  eq(4).  Similarly,  the  expectations  in  the  second  summation 
term  of  eq(8)  become 

E[yi(m)yJ(m)]  =  E[Vm]  r|)(0)  (13a) 

and 

E[y*(p)yj(p)]  =  E[v^)  [r;;®)]*.  (13b) 

Using  eqs(12)  and  (13)  in  (8), 
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1  K  K 

^rs  =  E  E  E[VmVp]  { lrj(0)|2  +  ru(m-p)[r]|(m-p)]*  } 

m=lp=l 

-  i^iiE[v^]E[Vp^lr;;(0)P 

^  m=lp=l 

1  K  K 

=  ^ E  E  {E[VmVp]  ru(m-p)[r]5(m-p)]*  +  D  lr”(0)|2}  (14) 

m=lp=l 

where 


D  =  E[vS,Vp]  -  E[v^]  E[vJ]  m,sp 
=  E[v<]  -  [E[v2]P  m  =  p 


(15a) 

(15b) 


Defining  q=m-p,  eq(14)  becomes 


K-l 


Vrs  =  ^  E  [K-lql]{E[VmVm.q]  ry(q)[r]5(q)]*  +  D  ^(0)|2}  m  =  1,2,. ..,K.  (16) 

^  q=-(K-l) 


Eq(16)  can  be  rewritten  as 


V,  =  ^{  EMr5(0)[r];(0)]*  +  {Elv"]  -  E^lv^] )  lrJ(0)P } 


1 

K 


K-l 


q=-(K-l) 
q^ 


1  - 


jql’ 

K 


{E[v^vLq]  r^(q)[r];(q)]*  +  D  lr;;(0)|2} 

m  =  1,2,...,K 


(17) 


where  the  first  term  is  the  q=0  (i.e.  the  m=p)  term.  Eq(17)  is  a  general  expression 
for  the  variance  of  the  sample  matrix  estimator  in  non-Gaussian  SIRP  disturbance. 
We  note  that  the  summation  term  involves  the  correlation  of  the  disturbance 
processes  over  the  range  cells. 


In  the  special  case  where  v  is  statistically  independent  over  range  cells, 
E[VmvLq]  =  E[Vm]  E[Vm J  m=l,2,...,K  (q^^O).  (18) 
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From  eq(15a),  D=0  for  this  case.  Further,  if  the  Gaussian  processes  Xj(k)  are 
uncorrelated  (and  therefore  statistically  independent)  over  range  cells, 

where  denotes  the  charmel  variance  associated  with  Xj(k)  and  will  be  assumed 
equal  on  all  channels  for  simplicity.  When  eq(19)  holds,  all  terms  in  the 
summation  of  eq(17)  are  zero,  so  that 

Vrs  =  H;{EMOx+  {ECv-tJ-EEv^J^llr^COiq.  (20) 

We  note  that  the  quantities  Ox  rjj(O)  relate  to  the  statistics  of  the 
processes  x(k)  which  are  not  directly  measurable.  However,  recognizing  that 
ay=E[v2]^ax  and  lifj(0)|2=E[v2]^lry(0)F,  we  can  easily  express  eq(20)  in  terms  of  the 
statistics  of  the  observed  data  processes  y(k)  with  the  result  that 

V„  =  i  +  {  E[v<]  -E[v2]^) }  lii;(0)|2 1 .  (21) 

We  note  that  the  variance  of  this  estimator  converges  to  zero  as  K 
approaches  infinity.  Thus,  the  sample  matrix  (using  spatial  averaging)  is  a 
consistent  estimator  for  non-Gaussian  SIRP's.  However,  this  quantity  is  now 
dependent  upon  the  second  and  fourth  moments  of  v. 

We  now  consider  two  special  cases  of  the  non-Gaussian  SIRP  models; 
namely,  the  Weibull  and  the  K-distribution.  For  K  distributed  processes  with 
shape  parameter  a  and  scale  parameter  b,  E[v4]=a(a-i-l)b'*  and  E[v2]=a/b2  (see 
Appendix  A).  We  let  E[y^]=l  without  loss  of  generality  so  that  E[v'^]=  1-i-l/a.  For 

XX  2 

simplicity,  consider  a  diagonal  element  term  [i.e.,  i=j  so  that  iij(0)=rii(Q)=ax].  In 
this  case  eq(20)  becomes 
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For  Gaussian  processes,  a->oo  so  that  eq(22)  reduces  to  the  well  known 

4 

result  a^/K.  However,  for  K-distributed  processes  with  high  tails  (a<l),  the  1/a 

term  causes  a  significant  increase  in  the  variance  of  the  estimator.  For  example, 
with  a=0.5,  is  five  times  larger  than  the  corresponding  Gaussian  case.  For 

WeibuU  processes  with  E[v2]=l  (see  Appendix  B), 


0<a<2 


For  the  case  of  WeibuU  processes  with  Rayleigh  amplitudes,  a=2  and  eq(23) 
also  reduces  to  a^/K.  For  WeibuU  processes  with  high  tails  (i.e.  low  a  values), 
increases  rapidly  with  decreasing  a.  FinaUy,  we  note  that  is  inversely  related 
to  the  beta  function  p(n,m)  with  m=n=2/a  where  p(n,m)=r(m)r(n)/r(m+n). 

2.2  Time-Averaged  Estimator 

We  now  consider  the  time-averaged  covariance  estimator;  i.e.  averaging 
performed  over  N  time  samples.  In  this  case, 

1  N-I5I 

=  (24) 

n=l 

where  6  denotes  temporal  lag.  Eq(8)  is  now  rewritten  as 

1  N-I5IN-I8I  *  ^ 

Vrt  =  E  E  E  {  yi(n)yj  (n-5)yi  (t)yj(t-5) } 

n=l  t=l 

1  N-I5IN-I5I  *  * 

-  ^  E  E  E[yi(n)yj  (n-6)]E[yi  (t)yj(t-5)]  (25) 

n=l  t=l 
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where  n  and  t  denote  the  i#’’  and  t*  time  samples,  respectively,  and  the  subscript  t 
on  V  denotes  the  time-averaged  case.  Again,  using  the  representation  theorem 
from  eq(4),  but  now  noting  that  v(k)  has  no  temporal  dependence,  the  expectations 
in  eq(25)  become 

E  {  yi(n)yj* (n-5)y * (t)yj(t-6) }  =  E[\^]  E  {  Xj(n)xj  (n-5)xi  (t)Xj(t-5) } .  (26a) 

E[yi(n)y*  (n-5)]  =  E[v2]E[Xi(n)xJ  (n-5)]  (26b) 

and 

E[y* (t)yj(t-5)]  =  E[v2]  E[Xi  (t)Xj(t-5)] .  (26c) 

Again,  for  zero  mean,  jointly  Gaussian  scalar  random  variables  Xj(*),  Xj(*), 
and  assuming  temporal  stationarity,  it  follows  that 

1  N-I5IN-I5I 

v„  =  S  EM  { lr^(5)P  +  rjl(n-t)[r];(n-t)]*  } 

n=l  t=l 

1  N-I5IN-I5I 
n=l  t=l 

=  R  ?  [l-^][EMru(v)[r|(v)]»]+ {E[v'']-E[v2]2)|r;;(5)P.  (27) 

v=-(N-l5l) 

where  v=  n-t  and  we  have  used  the  fact  that  l^j(5)I^  is  independent  of  v.  We 
emphasize  that  the  second  term  in  eq(27)  is  independent  of  N.  Thus,  the  variance 
of  the  time-averaged  estimator  does  not  converge  to  zero  for  increasing  N  which 
implies  that  this  estimator  lacks  consistency.  This  result  can  be  explained  by 
noting  the  non-ergodic  nature  of  the  SIRP  clutter  process  described  in  section  I. 

In  the  special  case  for  i=j  and  5=0,  expresses  the  variance  Var(6^2)  of  the 
estimated  power  on  a  given  channel  as  [8] 

Var(ft2)=^£  [l-^]EMIrjl(v)P+{E[v'']-E[v2]^)o"  (28) 

v=-N 
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Furthermore,  for  white  noise  processes,  only  the  v=0  term  in  the  summation 

X  2 

remains,  so  that  with  rji(0)=Ox 

Var(^2)  =  1  E[v4]ax  +  {E[v^]  -  E[v2]^}  (29) 

which  compares  to  the  result  presented  in  [3].  Again,  we  note  the  inconsistency  of 
the  estimator.  Using  the  previously  defined  expressions  of  E[v^]  and  E[v‘^],  for  K- 
distributed  processes 

Var(&2)  =  ^  [l  +  ^]  oj  +  ^  (30) 


while  for  Weibull  processes 


(31) 


For  the  Gaussian  case,  <x=oo  and  oc?=2  in  eqs(30)  and  (31),  respectively,  and 
Var(6^2)  reduces  to  the  known  result  a^/N. 
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3.0  Simulated  Results 


Validation  of  eqs(22)  and  (23)  was  obtained  via  Monte-Carlo  using  1000 

4 

statistically  independent  realizations.  The  normalized  term  KV/Ox  is  plotted  as  a 

function  of  the  shape  parameter  a  for  the  K-distributed  noise  and  the  Weibull  noise 
in  Figures  1  and  2,  respectively.  These  curves  approach  unity  (the  Gaussian  case) 
as  a-»oo  in  Figure  1  for  K-distributed  noise  and  at  a=2  (Rayleigh  amplitudes)  in 
Figure  2  for  Weibull  noise.  These  curves  show  the  dramatic  increase  in  the 

4 

variance  as  a  function  of  the  shape  parameter.  In  Figure  3,  Vj/a^  is  plotted  as  a 
function  of  K,  the  number  of  independent  range  cell  samples,  for  the  Gaussian  and 

4 

K-distributed  (a=0.1  and  0.5)  cases.  These  curves  show  that  V/ax  decreases 

linearly  with  increasing  K.  Thus,  the  estimator  is  consistent  when  averaging  over 

4 

range  cells.  This  is  unlike  the  time-averaged  covariance  estimator  in  which  V/Ox 

converges  to  a  bias  level  with  increasing  number  of  time  samples  due  to  the  non- 
ergodic  nature  of  the  SIRP's.  Finally,  they  reveal  the  number  of  samples  K 
required  to  reduce  the  variance  of  the  estimator  for  the  non-Gaussian  processes  to 
that  of  the  Gaussian  case. 

For  the  time-averaged  estimator,  we  plot  the  normalized  variance 
versus  N  for  eqs(30)  and  (31)  in  Figures  4  and  5,  respectively.  These  figures  reveal 
the  significance  of  the  inconsistent  time-averaged  covariance  estimator  in  non- 
Gaussian  SIRP's. 
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Normalized  Variance 


Shape  Parameter  for  K-Distribution 

Figure  1  Normalized  variance  K(Vj/a5  for  a  diagonal  element  of  the  covariance 
matrix  estimator  versus  the  shape  parameter  for  K-distributed  noise. 
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Normalized  Variance 


0  1  2 

Shape  Parameter  for  Weibull  Distribution 


Figure  2  Normalized  variance  K(Vj/a5  for  a  diagonal  element  of  the  covariance 

matrix  estimator  versus  the  shape  parameter  for  Weibull  distributed 
noise. 
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Normalized  Variance 


Number  of  Independent  Samples  (K) 


Figure  3  Normalized  variance  (V/a^  versus  the  number  of  independent  data 

samples  for  a.)  O  Gaussian  noise,  b.)  ^  K-distribution  (a=0.5), 

c.)  •  K-distribution  (a=0.1).  Simulated  experimental  values  shown  by 

symbols. 
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Normalized  Variance  of  Time-Averaged  Estimator 


Figure  4  Normalized  variance  for  the  time-averaged  correlation  matrix  estimator 
element  for  K-distributed  processes  versus  the  time  sample  window  size 
N  with  a  as  a  parameter. 
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Figure  5  Normalized  variance  for  the  time-averaged  correlation  matrix  estimator 
element  for  Weibull  distributed  processes  versus  the  time  sample  window 
size  N  with  a  as  a  parameter. 


4.0  Summary 


This  paper  presents  analytical  expressions  for  the  variance  of  the  the 
estimated  covariance  matrix  elements  when  the  observation  data  processes  are  non- 
Gaussian  spherically  invariant  random  processes  (SIRP's).  Both  the  spatially- 
averaged  and  time-averaged  are  considered.  These  expressions  show  the 
significant  increase  in  the  variance  of  the  estimates  as  compared  to  the  Gaussian 
case.  For  spatial  averaging,  they  indicate  the  increase  in  the  number  of  samples 
over  which  averaging  must  be  performed  to  acheive  the  same  level  of  performance 
as  the  Gaussian  case.  For  the  time  averaging,  the  non-ergodic  nature  of  the  SIRP 
clutter  precludes  convergence  of  the  estimator  variance.  Validation  is  acheived 
through  Monte-Carlo  simulation. 
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APPENDIX  A 


DERIVATION  OF  THE  EXPRESSION  FOR  EEv^]  USING 
K-DISTRIBUTED  PROCESSES 

In  this  Appendix,  an  analytic  expression  for  E[v^]  is  derived  in  terms  of  the 
shape  and  scale  parameters  of  ^(v)  where  fyCv)  is  a  Gamma  distribution  and  the 

random  variable  v  is  real  and  non-negative.  In  this  case,  the  observation  data 
vector  process  y(n,k)=v(k)x(n,k)  described  by  eq(4)  has  a  K-distributioh.  For  the 
Gamma  distributed  fyCv)  with  non-negative  v, 

fv(v)  =  (bv)2a-i  expC-b^v^) 

so  that 

oo 

E[v'^]  =  J  (bvp^'V'^  exp(-b^v2)dv 

i  (^)  0 

oo 

=  f  (b)2a-iv2“+3  exp(-b2v2)dv.  (A.l) 

r(a)  (/ 

Consider  v=VWb  so  that  dv/dw=l/2b'\/^  Now, 


fw(w)  = 


w”'^exp(-w) 

r(a) 


and 

E[v4]  =  E[w2]/b4. 

However, 

OO 

_r  f  w“~^exp(-w)  ,  r(a+2) 

E[w]-J  dw-  -a(a+l). 

Using  (A.4)  in  (A.3) 


(A.2) 


(A.3) 


(A.4) 


E[v'^]  =  E[w2]/b^  = 


(A.5) 
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APPENDIX  B 


DERIVATION  OF  THE  EXPRESSION  FOR  E[s^]  USING 
WEIBULL  DISTRIBUTED  PROCESSES 

In  this  Appendix,  an  analytic  expression  for  E[v^]  is  derived  in  terms  of  the 
shape  and  scale  parameters  for  processes  with  a  Weibull  distributed  amplitude. 
The  Weibull  probability  density  function  (pdO  of  the  amplitude  r  denoted  as  f(r) 
with  shape  parameter  a  and  scale  parameter  b  is  expressed  as 

f(r)  =  ba  r“'^exp(-br“)  r^,  0<a<2.  (B.l) 

The  conditional  pdf  f(rls)  is 

f(rlv)  =  ^  exp(-rV2v2)  r^.  (B.2) 

Using(B.l),  the  k*  moment  of  r  is 

oo 

E[r*^]  =  J  r‘^f(r)dr 
0 

oo 

=  ba  J  i^'^’^exp(-br“)dr 

=  (E|'“r(l+k/a)  (B.3) 

while  the  moment  of  r  conditioned  on  v  is  obtained  using  (B.2),  so  that 

OO 

E[r*^lv]  =  J  i^f(rlv)dr 
0 

Let  w=r2/2v^,  so  that  dw=  rdr/v^  and  r='\/2vw^^.  Eq(B.4)  now  becomes 


f  ";^exp(-r2/2v2)dr.  (B.4) 
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oo 


E[r^lv]  =  2*^  J  exp(-w)dw 

=  2*^v*^r(l+k/2).  (B.5) 


Now, 

OO 

E[r*^]  =  J  E[r*^lv]  f(v)dv 

OO 

=  2‘^r(l+k/2)  J  v‘'f(v)dv 

=  2k/2  r(l+k/2)  E[v*^] 

Using  (B.3)  in  (B.6), 

fsT"  r(i+k/a) 

E[v'‘]=  2kPr(l+k/2) 

Therefore, 


ForE[v2]=l, 


a  lo/^ 

r(2/a). 


(B.6) 


(B.7) 


(B.8) 


(B.9) 


so  that 


E^-i]  = 


I'g'i  r(4/a) 

y^h\2/ay 
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